function [R_S_GHQ,R_S_GHQ_RN,weight] = CreateNetReturns_GHQ(lambda, mu_d_RN, MRP, sigma_d, mu_j, sigma_j)

% Generate a Gauss-Hermite Quadrature of stock returns based on Merton Jump Diffusion (MJD) model.  
    % Output (stored): GHQ nodes for net return of market index:
        % 'R_S_GHQ'     (n-by-1 vector of net returns under real-world measure)
        % 'R_S_GHQ_RN'  (n-by-1 vector of net returns under risk-neutral measure)
        % 'weight'      (n-by-1 vector of GHQ weights)

mu_d = mu_d_RN + MRP;
tau = 1;            % length of return period (in years)
delta_div = 0;      % annual dividend rate


%% Gauss-Hermite Quadrature

N_P = max(10,round(5*lambda,0))*tau;    % max. number of jumps in a given period.
p_j_k = zeros(N_P+1,1); % Poisson prob. of 'k' jumps
for k = 0:N_P
    p_j_k(k+1) = exp(-lambda*tau)*(lambda*tau)^k/prod(1:k);
end
N_G = 500;      % number of Hermite nodes
[xx,ww] = GaussHermite(N_G); % xx: column vector of Hermite nodes, ww: column vector with corresponding weights.
xx = sqrt(2)*xx;    % adjustment to reflect N(0,1) distribution
ww = ww/sqrt(pi);   % adjustment to reflect N(0,1) distribution

LR = zeros((N_P+1)*N_G,1);      % GHQ sample log-returns of index over crediting term
LR_RN = zeros((N_P+1)*N_G,1);	% GHQ sample log-returns of index over crediting term, under risk-neutral measure
weight = zeros((N_P+1)*N_G,1);  % GHQ sample weights
for k = 0:N_P
    LR( (k*N_G+1):((k+1)*N_G) ) = ( mu_d - delta_div - sigma_d^2/2 )*tau + k * mu_j + sqrt( sigma_d^2 * tau + k * sigma_j^2 )*xx;
    LR_RN( (k*N_G+1):((k+1)*N_G) ) = ( mu_d_RN - delta_div - sigma_d^2/2 )*tau + k * mu_j + sqrt( sigma_d^2 * tau + k * sigma_j^2 )*xx;
    weight( (k*N_G+1):((k+1)*N_G) ) = p_j_k(k+1) * ww;
end

% pruning:
pruning_threshold = 1e-7;
LR = LR( (weight>pruning_threshold) & (weight<Inf) );
LR_RN = LR_RN( (weight>pruning_threshold) & (weight<Inf) );
weight = weight( (weight>pruning_threshold) & (weight<Inf) );

R_S_GHQ = exp(LR) - 1;  % representative sample of net index returns over crediting term.
R_S_GHQ_RN = exp(LR_RN) - 1;    % representative sample of (risk-neutral) net index returns.

return